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Abstract 

r-H 

This paper puts forth a new large eddy simulation closure modeling strategy for two-dimensional turbu- 
lent geophysical flows. This closure modeling approach utilizes approximate deconvolution, which is based 

model. The new approximate deconvolution model is tested in the numerical simulation of the wind-driven 
circulation in a shallow ocean basin, a standard prototype of more realistic ocean dynamics. The model 
employs the barotropic vorticity equation driven by a symmetric double-gyre wind forcing, which yields a 
four-gyre circulation in the time mean. The approximate deconvolution model yields the correct four-gyre 



circulation structure predicted by a direct numerical simulation, on a coarser mesh but at a fraction of the 
computational cost. This first step in the numerical assessment of the new model shows that approximate 
deconvolution could represent a viable tool for under-resolved computations in the large eddy simulation of 
more realistic turbulent geophysical flows. 

ft 
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1. Introduction 

o 

r-H 

As a first approximation, the mean near-surface currents in the oceans are driven by the mean effects 

• • 

^ of winds. Wind-driven flows of mid-latitude ocean basins have been studied by modelers using idealized 

^ single and double-gyre wind forcing. This type of double-gyre circulation characterizes all mid-latitude 

ocean basins, including those in the North and South Atlantic, as well as the North and South Pacific. 
The barotropic vorticity equation (BVE) represents one of the most commonly used mathematical models 



for this type of geostrophic flows with various dissipative and forcing terms (Majda and Wang 2006). For 



more details on the physical mechanism and formulations of the BVE, the reader is referred to Holland 
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and Rhines|(Tl980); Munk and Wunschl (119821); IGriffa and Salmonl (119891); ICumminsI (119921); IGreatbatch and 



Nadiga (2000); Nadiga and Margolin (2001); Fox-Kemper (2005); Cushman-Roisin (1994). 



The main assumptions that go into the BVE model are hydrostatic balance, the /3-plane approximation, 
geostrophic balance, and horizontal eddy viscosity parametrization. Despite the fact that the BVE models 
are a simplified version of the full-fledged equations of geophysical flows, their numerical simulation is still 
computationally challenging when long-time integration is required, as is the case in climate modeling. 
Thus, to further reduce the computational cost of the BVE, large eddy simulation (LES) appears to be a 
natural approach. In LES only the large spatial structures are approximated, whereas the small scales are 
modeled. This allows for much coarser spatial meshes and thus a computational cost that is significantly 
lower than that of a direct numerical simulation (DNS). To achieve the same order of physical accuracy as 



DNS, however, LES needs to correctly treat the closure problem (Meneveau and Katz, 2000 Sagaut, 2006 



Berselli et al. , 2006): the effect of the small scales on the large ones needs to be modeled. This closure 



problem represents a significant challenge for quasi- two-dimensional turbulent flows, such as those in the 



ocean and atmosphere ( McWilliams , 1989 Danilov and Gurarie, 2000 Smith et al. , 2002 



Lesieur 



2008 



Salmon, 1998 Majda 2003 Holton 2004 McWilliams, 2006 Vallis, 2006). 



Two-dimensional turbulence is a fundamental topic for understanding of geophysical flows and behaves 
in a profoundly different way from the three-dimensional turbulence due to different energy cascade behavior 



(Kraichnan, 1967 Batchelor 1969 Leith 1971), which is described in the Kraichnan-Batchelor-Leith two- 
dimensional turbulence theory. In three-dimensional turbulence, energy is transferred forward, from large 
scales to smaller scales, via vortex stretching. In two dimensions that mechanism is absent, and under most 
forcing and dissipation conditions energy is transferred from smaller scales to larger scales, largely because of 
the potential enstrophy, a quadratic invariant defined as the integral of the square of the potential vorticity. 
The physical mechanism behind the enstrophy cascade is the stretching of small-scale vorticity gradients 



by the strain arising from larger-scale vortices (Chen et al. , 2003). This results in energy being trapped 
at large scales in forced-dissipative two-dimensional flows. Additionally, enstrophy is transferred to the 



smaller scales from larger scales, and is finally dissipated at large wavenumbers. Danilov and Gurarie (2000) 



and Tabeling (2002) reviewed both theoretical and experimental two-dimensional turbulence studies and 
provided extensive insights into the applicability of two-dimensional turbulence theory to geophysical flows. 



We also refer the reader to the energy-enstrophy conservation arguments provided by |Vallis (2006). 

Modeling the ocean and atmosphere inspired the first LES models, but most of the subsequent devel- 



opment of LES has taken place in the engineering community (Sagaut, 2006 Berselli et al. 2006). The 
majority of LES models have been developed for three-dimensional turbulent flows, such as those encoun- 
tered in engineering applications. These LES models fundamentally rely on the concept of the forward 
energy cascade and so their extension to geophysical flows is beset with difficulties. The effective viscosity 
values in oceanic models are much greater than the molecular viscosity of seawater, hence a uniform eddy 

2 



viscosity coefficient is generally used to parameterize the unresolved, subfilter-scale effects in most oceanic 



models ( McWilliams , 2006 Vallis, 2006). LES models specifically developed for two-dimensional turbulent 



flows, such as those in the ocean and atmosphere, are relatively scarce (Fox-Kemper and Menemenlis , 2008 



A wad et al. 



2009 



Ozgokmen et al. 



2009 



Chen et al. 



2010), at least when compared to the plethora of 



LES models developed for three-dimensional turbulent flows. |Holm and Nadi ga (2003) combine the uniform 
eddy viscosity parametrization with the alpha regularization LES approach to capture the under-resolved 
flow where the grid length becomes greater than the specified Munk scale of the problem. In that work the 
structural alpha parameterization was tested on the BVE in an ocean basin with double-gyre wind forcing, 
which displays a four-gyre mean ocean circulation pattern. It was found that the alpha models provide a 
promising approach to LES closure modeling of the barotropic ocean circulation by predicting the correct 
four-gyre circulation structure for under-resolved flows. 

This paper puts forth a new LES closure modeling strategy for two-dimensional turbulent geophysical 
flows. The new closure modeling approach utilizes approximate deconvolution (AD), which is particularly 
appealing for geophysical flows because of no additional phenomenological approximations to the BVE. The 
AD approach can achieve high accuracy by employing repeated filtering, which is computationally efficient 
and easy to implement. The AD method has been used successfully in LES of three-dimensional turbulent 
engineering flows ( |Stolz and Adam s, 1999; Stolz et al. , 2001a|b 2004). We emphasize, however, that to 
the best of our knowledge, this is the first time that the AD methodology is used in LES of large scale 
geophysical flows, such as the barotropic ocean circulation flow we consider in this paper. To assess the 
new AD closure modeling approach, we test it on the same two-dimensional barotropic flow problem as that 



employed in Nadiga and Margolin (2001) and in Holm and Nadiga (2003). 



The rest of the paper is organized as follows: The BVE, the mathematical model used in this report, is 
presented in Section [2] Section [3] presents the AD methodology and introduces the new closure model. The 
numerical methods used in our simulations are briefly discussed in Section [4j The results for the new AD 
model are presented in Section [5] Finally, the conclusions are summarized in Section [6] 

2. Barotropic Vorticity Equation 



In this section, we present the BVE, the mathematical model used in the numerical investigation of the 
new AD model. The BVE is one of the most used mathematical models for geostrophic flows with various 



dissipative and forcing terms (Majda and Wang, 2006). Studies of wind-driven circulation using an idealized 



double-gyre wind forcing have played an important role in understanding various aspects of ocean dynamics, 



including the role of mesoscale eddies and their effect on mean circulation. Following Greatbatch and Nadiga 



(2000), we briefly describe the BVE. For more details on the physical mechanism and various formulations 



utilized, the reader is referred to Holland and Rhines (1980); Munk and Wunsch (1982); Griffa and Salmon 
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( 1989[); |Cummins| fl!992|); |Greatbatch and Nadiga] (|2000|); [Nadiga and Margolin| (|200T|) ; | Fox- Kemp er| d2QQ5[). 



The governing equations for two-dimensional incompressible barotropic flows can be written in dimen- 
sionless form of the potential vorticity formulation in the beta plane as the BVE: 

dq 



dt 



J = D + F. 



In Eq. ([I]), q is the potential vorticity, defined as 

q = Rouo + y, 



(1) 



(2) 



(3) 



where uo is the vorticity and Ro is the Rossby number. 

The nonlinear convection term in Eq. ([!]), called the Jacobian, is defined as 

dq dip dq 
dy dx dx dy ' 

where ip is the streamfunction. The kinematic relationship between the vorticity and the streamfunction 
yields the following Poisson equation: 

(4) 



d 2 ^ d 2 ^ _ 
dx 2 dy 2 

The viscous dissipation in Eq. has the form 



5m\ 3 ( d 2 u d 2 u 



^ 1 L J \ dx 2 ' dy 2 J ' 
where 5m is the Munk scale and L is the basin dimension. The double-gyre wind forcing is given by 

F = F sm(iry), 



(5) 



(6) 



where Fq = 1 due to the Sverdrup velocity scale used for nondimensionalization (Greatbatch and Nadiga 



2000). In this nondimensionalization, the velocity scale is 

7T T 



V 



(7) 



where ro is the maximum amplitude of double-gyre wind stress, p is the mean density, H is the mean depth 
of the basin, and j3 is the gradient of the Coriolis parameter at the basin center (y = 0). The dimensionless 
variables are then defined as 

t 



x 

x = -, y 



y 



t 



o q 



(8) 



V ° V L/V 

where the tilde denotes the corresponding dimensional variables. In dimensionless form, there are only two 
physical parameters, the Rhines scale and the Munk scale, which are related to the physical parameters in 
the following way: 

A, / V X 1 / 2 K., f u \!/3 

(9) 
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where v is the uniform eddy viscosity coefficient. The physical parameters in the BVE ([T]), the Rhines scale 
5 1 and the Munk scale 5m, are related to the Reynolds and Rossby numbers through the following formulas: 



L 



(Ro) 



1/2 



(10) 



°M 

L 



(Re _1 Ro) 



1/3 



(11) 



where Re is the Reynolds number based on the basin dimension, L. We note that some authors use a 
boundary layer Reynolds number, which is written as 

Si _ 5 3 j 



Re# = Re 



L 



5 3 ' 

°M 



(12) 



where Re# ~ O(10) — O(10 3 ) for oceanic flows (Fox-Kemper, 2005). Finally, in order to completely specify 
the mathematical model, boundary and initial conditions need to be prescribed. In many theoretical stud- 
ies of large scale ocean circulation, slip or no-slip boundary conditions are used. Following these studies 



(e.g., Greatbatch and Nadiga, 2000 Nadiga and Margolin 2001 Holm and Nadiga 2003 Cummins, 1992 



Ozgokmen and Chassignet 



1998 



Munk 



1950 



Bryan 1963), we use slip boundary conditions for the velocity, 



which translate into homogenous Dirichlet boundary conditions for the vorticity: uj\q = 0. The imperme- 
ability boundary condition is imposed as = 0. For the initial condition, we start our computations from 
a quiescent state (q = 0) and integrate Eq. ([I]) until a statistically steady state is obtained in which the 
wind forcing, dissipation, and Jacobian balance each other. 



3. Approximate Deconvolution Model 

The AD approach aims to obtain accurate and stable approximations of the original, unfiltered flow 



variables when approximations of the filtered variables are available (Stolz and Adams, 1999 Germano 



2009). The AD methodology was developed in the image processing community, and has been successfully 
adapted to the closure problem in turbulence modeling for engineering flows (Stolz et al. , 2001a|b Adams and 



Stolz 2002 Stolz et al. , 2004). We emphasize that this approach is purely mathematical, with no additional 



phenomenological arguments being used. This is particularly appealing for LES of geophysical flows in which 
different energy transfer characteristics are displayed than those in three-dimensional turbulent flows, for 
which successful phenomenological modeling has been done. Next, we present the mathematical derivation 
of the new AD model for the BVE 0. 

To derive the equations for the filtered flow variables, the BVE ([I]) is first filtered with a rapidly decaying 
spatial filter (to be specified later). Thus, using a bar to denote the filtered quantities, the filtered BVE 
reads: 

dq 
di 

5 



■J(q,$) = D + F + S, 



(13) 



where S is the subfilter-scale term, given by 

S = 



-J(q,iP) + J(q,iP). (14) 
It is precisely at this point in the LES model derivation that the celebrated closure problem must be 



addressed. In order to close the filtered BVE (13), the subfilter-scale term S in Eq. (14) needs to be 
modeled in terms of the filtered flow variables, q and t/S. 

This paper proposes a new LES closure modeling approach for two-dimensional turbulent geophysical 
flows. The AD approach can achieve high accuracy, is computationally efficient, and is easy to implement. 
Although the AD methodology has already been successfully used in LES of three-dimensional turbulent 
engineering flows (Stolz and Adams 1999 Stolz et al. , 2001a|b 2004), this is the first time that it is used in 
LES of large scale geophysical flows, such as the barotropic ocean circulation flow considered in this paper. 

The goal in AD is to use repeated filtering in order to obtain approximations of the unfiltered flow 
variables when approximations of the filtered flow variables are available. These approximations of the 
unfiltered flow variables are then used in the subfilter-scale terms to close the LES system. To derive the 
new AD model, we start by denoting by G the spatial filtering operator: Gq = q. Since G = I — (I — G), an 
inverse to G can be written formally as the non- convergent Neumann series: 



g- 1 -^Tii-Gy- 1 . 



(15) 



Truncating the series gives the van Cittert approximate deconvolution operator, Qat. We truncate the series 
at N and obtain Qn as an approximation of G~ x : 

N 

Q N = J2(i-Gy-\ (i6) 

i=l 

where I is the identity operator. The approximations Qn are not convergent as N goes to infinity, but rather 



are asymptotic as the filter radius, (5, approaches zero (Berselli et al. 2006). An approximate deconvolution 
of q can now be obtained as follows: 

q * q* = QnQ- (IT) 



For higher values of N we get increasingly more accurate approximations of q: 



Qi 

Q2 

Q 3 
Qs, 



I 

21 -G 

3/ - 3G + G 2 

4I-QG + 4G 2 - G 3 

51 - 10G + 10G 2 - 5G 3 + G 4 



(18) 
(19) 
(20) 
(21) 
(22) 
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Following the same approach as that used in Dunca and Epshteyn ( 2006| ), one can prove that these models 
are highly accurate (0(5 2N + 2 ) modeling consistency error) and stable. We choose N = 5 and find an AD 
approximation of the variable q as 



q « q* = 5q - lOq + 10 q - bq + q 



and, similarly, an AD approximation of the variable ip\ 



i/j « ^* = 5^ - 10^ + 10^ - 5^ + i/i- 



Using (23) and (24), we can now approximate the nonlinear Jacobian: 



J(q^) « J(g*,^*). 



(23) 



(24) 



(25) 



Finally, using the approximation Eq. (25) in Eq. (14), closes the filtered BVE (13) and yields the new AD 
model: 



dq 
dt 



■J(q^) = D + F + S* 



where S* is the subfilter-scale term, given by 



S* = -J(q*,ip*) + J(q,ip)- 



(26) 



(27) 



To completely specify the new AD model (26)-(27), we need to choose a computationally efficient filtering 



operator. Following Stolz and Adams (1999), we use the following second-order accurate filtering operator: 

fi-1 + fi+1 



i-l 



fi + Ctfi+l 



fi 



(28) 



where the subscript i is the spatial index in the ^-direction. This results in a tridiagonal system of equations 



for each fixed value of y. A generalization of Eq. p8j), the spectral- type high-order compact filter introduced 
in f 

M 



Visbal and Gaitonde ( 2002 ) , is given by 

a/i-i + f% - 



■ Oifi+i — ~^(fi+n + fi-n), 



(29) 



n=0 



where / is the computed value, and / is the corresponding once filtered value. This filter provides 2Mth- 



order accurate filtering on a 2M + 1 point stencil (Visbal and Gaitonde, 2002). The free parameter, a, which 



is in the range < \a\ < 0.5, determines the filtering properties, with high values of a yielding less dissipative 
results. The coefficients for different order filters are given in Table [l] In our numerical tests, we found 
that for moderately fine meshes, the range of 0.25 < a < 0.5 is appropriate, but for coarser meshes, lower 
values (e.g., a = 0.1) should be used to eliminate spurious oscillations. Most of our numerical simulations 



use the free parameter a = 0.25 (Stolz and Adams 1999). Since we used a second-order accurate space 



discretization in this study, the filtered simulations are conducted either by using the second-order accurate 
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2nd-order 4th-order 6th-order 8th-order 
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2 
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16 
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14a 
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8 


32 ~ 


32 


as 










1 

32 


OL 

16 


1 

16 


a 
8 


(I4 















1 

128 


1 OL_ 

^ 64 



Table 1: Coefficients of the compact niters of various orders given by Eq. \29\ 



filter given in Eq. (28) or the fourth-order accurate compact filter given in Eq. (29) with different parameter 
values for a. 

The existence and uniqueness of strong solutions of the BVE model using the AD closure was investigated 



by Stanculescu (2008). According to this work, the transfer function of the filter used in the AD closure 



should be positive definite. The transfer function of the tridiagonal filter given in Eq. (28) can be written 

as 

I 1 -I- rnsf f,A 

(30) 



TM = (i + a) r 1 + C ° SM 



■ 2a cos(cj) 

The transfer function of the generalized filter given in ([29]) can be written as 



2^ n =0 a n COS ( n 1 



T(u>) 



(31) 



1 + 2a cos(co>) 

The transfer functions corresponding to the second- and fourth-order filters are plotted in Fig. ([!]) and are 
positive definite in the interval of < \a\ < 0.5, where a is a free parameter of the filter. By substituting 
these transfer functions into the deconvolution operator, it can be easily seen that the transfer functions of 
the deconvolution operators also become positive definite. With these two conditions verified (the positive 
definiteness of the transfer functions of the filters and deconvolution operators), the theoretical results in 



Stanculescu (2008) provide solid mathematical support for our numerical investigations. 



4. Numerical Methods 

In this section, we provide a brief description of the numerical methods employed in the investigation of 
the new AD model. For the time discretization, we employ an optimal third-order total variation diminishing 



Runge-Kutta (TVDRK3) scheme (Gottlieb and Shu, 1998). For clarity of notation, we rewrite the new AD 
model in the following form 




9 



where R denotes the discrete spatial derivative operator, including the nonlinear convective term, the linear 
diffusive term, the forcing term, and the subfilter-scale term. The TVDRK3 scheme then becomes 



q 



(i) 



=7(2) 



q n+1 



• AtR M 



r(2) 



l -AtR( 2 



(33) 



3" 3" 

where At is the adaptive time step size which can be computed at the end of each time step by: 

min(Ax, Ay) 



At 



max 



(if uffi)' 



(34) 



where c is known as Courant-Friedrichs-Lewy (CFL) number. We set c = 1.0 for all the simulations reported 
here. 

For the spatial discretization of linear operators, such as the viscous dissipation term and the Poisson 
equation, we use a second-order finite difference method. For the nonlinear convection term, we utilize 



Arakawa's conservative scheme (Arakawa 1966). Finally, a direct Poisson solver based on the fast sine 
transform is used to solve the kinematic relationship between the vorticity and the streamfunction in Eq. Q. 
Specifically, we use a fast sine transform in one direction and a tridiagonal system solver in the other direction 



(Moin 2001). 



For completeness, we summarize the solution algorithm for one time step. We assume that the numerical 
approximation for time level n is known, and we seek the numerical approximation for time level n + 1, 



after the time step At. We write the algorithm for the new AD model (26)-(27), but we emphasize that it 
is the same algorithm, with some obvious modifications, that is used for the high resolution well resolved 
simulation, which is referred to here as DNS. We emphasize that the term DNS in this study is not meant 
to indicate that a fully detailed solution is being computed, but instead refers to resolving the simulation 
down to the Munk scale via the specified lateral eddy viscosity parameterization. The solution algorithm 
has the following form 

1. Compute tp n using the Poisson solver for the following equation: 

Q2^n Q2^n X 



dx 2 



dy 2 



Ro 



(q n - y) • 



(35) 



2. Compute q n * and ij) n * using q n and the filtering procedure in Eq. (29), and the AD method in 



Eq. (23) and Eq. (24). 



3. Compute the nonlinear Jacobian terms J(q n *,tjj n *) and J(q n ^ n ) using the Arakawa scheme. 



4. Compute the subfilter-scale term S* using J(^ n *, ^ n *), J(g n , ip n ), and Eq. (27). 
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5. Compute the viscous dissipation term D using a central difference scheme and calculate the right- 



hand-side term R in Eq. (32). 



6. Apply the Runge-Kutta algorithm by repeating steps 1-5 for each stage of Runge-Kutta scheme. 



7. Adjust the time step according to the CFL number given in Eq. (34). 

8. Go to step 1 to begin the next time step. 



5. Numerical Results 



The main goal of this section is to test the new AD model (26)-(27) in the numerical simulation of the 



BVE ([!]). Before testing the new model, however, we first validate the code in Section 5.1 The new AD 
model is then tested in the numerical simulation of a four-gyre barotropic model with double-gyre wind 



forcing in Section 5.2 



5.1. Method of manufactured solutions: Code validation 

To validate the code, we consider the following exact solution, which corresponds to a steady state 
Taylor- Green vortex flow: 

LJ 

2 7r 2 sin(7nr) sin(7n/), 
sin(7rx) sin(7n/). 



2 7r 2 sin(7rx) sin(7n/) + y, 



(36) 
(37) 
(38) 



The forcing function is calculated by plugging Eqs. (36)-(38) into the BVE 0: 



F(x,y,t) 



-7r cos(7rx) sm(iry) + ( ) 47r 4 sin(7nr) sm(ny). 



(39) 



We use q(x, y, 0) =0 as the initial condition. The spatial resolution is 64 x 128 and the time step At is 
computed adaptively as given in Eq. ( 34 ) . We utilize the numerical algorithm in Section [4] to integrate the 
BVE until the steady state solution is obtained. Our numerical results closely match the exact solution. To 
investigate the individual effects of the nonlinear Jacobian, dissipation, forcing and subfilter-scale terms, we 
compute the time series of following integral quantities: 



Qj(t) = \ Jl (J(q,t)) 2 dxdy, 



(40) 



?o(t) = \ II D 2 dxdy. 



?f(() = g jj F 2 dxdy, 
11 



(41) 
(42) 




Figure 2: Code validation test case. Time history of the total energy, individual dissipation, Jacobian, subfilter-scale and forcing 
terms: (a) 8j/L = 0.04 and 5m /L = 0.02; and (b) Sj/L = 0.06 and 5 m /L = 0.02. Note that the numerical approximation 
converges to the correct value after a short transient interval. 



h(t) = l 1 1 (S*) 2 dxdy. (43) 



We also plot with respect to time the total energy: 

ill 



m-i it + (44) 



2 J J \dy J \dx 



Since the exact steady state solution is given in Eqs. (36)- (38), we can calculate the exact values of the 
above integral quantities. 

In Fig. ([2| we plot the time evolution of the total energy, and contributions of the individual dissipation, 
subfilter-scale, nonlinear Jacobian, and forcing terms. We plot these time series for two parameter sets: 
5i/L = 0.04 and 5 M /L = 0.02 in Fig. $2^), and <S 7 /L = 0.06 and 5 M /L = 0.02 in Fig. ([§)). For both 
parameter sets, the time evolution of the above integral quantities follows the same pattern: after a short 
transient interval, they converge to the correct exact values. Thus, we conclude that the numerical algorithm 
produces accurate results for these validation test cases. 

5.2. Four- gyre problem with double- gyre forcing 



The main goal of this section is to test the new AD model Eqs. (26)- (27) in the numerical simulation of 
the wind-driven circulation in a shallow ocean basin, a standard prototype of more realistic ocean dynamics. 
The model employs the BVE driven by a symmetric double-gyre wind forcing, which yields a four-gyre 



circulation in the time mean. This test problem has been used in numerous studies (e.g., Cummins 1992 
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Greatbatch and Nadiga] |2QQQ| |Nadiga and Margolin) |2001| |Holm and Nadiga[ |2QQ3| |Fox-Kemper] |2005[). 



This problem represents an ideal test for the new AD model. Indeed, as showed in |Greatbatch ancLNa diga 



( 2000| ), although a double gyre wind forcing is used, the long time average yields a four gyre pattern, which 
is challenging to capture on coarse spatial resolutions. Thus, we will investigate numerically whether the 
new AD model can reproduce the four gyre time average on a coarse mesh. 

The mathematical model used in the four gyre problem is the BVE ([!]). The double-gyre forcing is 
given in Eq. (J6|. Following |Holm and Nadi ga (2003), we utilize two different parameter sets, corresponding 
to two physical oceanic settings: Experiment (i) with a Rhines scale of Sj/L = 0.04 and a Munk scale of 
5m /L = 0.02, which corresponds to a Reynolds number of Re = 200 (or a boundary layer based Reynolds 
number of Rcb = 8) and a Rossby number of Ro = 0.0016; and Experiment (ii) with a Rhines scale of 
Si/L = 0.06 and a Munk scale of 5 m /L — 0.02, which corresponds to a Reynolds number of Re = 450 
(or a boundary layer based Reynolds number of ReB = 27) and a Rossby number of Ro = 0.0036. Since 
we set the Munk scale to 5m /L = 0.02 in our study, the uniform eddy viscosity coefficient embedded in 
the model can be calculated from Eq. For example, if we set the mid-latitude ocean basin length to 
L = 2000 km and the gradient of the Coriolis parameter to f) = 1.75 x 10 -11 m _1 s _1 , then our model uses 
v = 1120 m 2 s _1 as its eddy viscosity parametrization. All numerical experiments conducted here are solved 
for a maximum dimensionless time of T max = 100. This value corresponds to the dimensional times of 
56.6 and 25.15 years for Experiment (i) and Experiment (ii), respectively, which are long enough to capture 
statistically steady states. All computations were carried out using the gfortran compiler on a Linux cluster 
system (3.4 GHz/node). 

To assess the new AD model, we employ the standard LES methodology: we first run a DNS computation 
on a fine mesh. We then run an under-resolved numerical simulation on a much coarser mesh (denoted in 
what follows as BVE coarse ), which does not employ any subfilter-scale model. We expect BVE coarse to 
produce inaccurate results. Finally, we employ the new AD model on the same coarse mesh utilized in 
BVE coarse . The AD model should yield results that are significantly better than those obtained with 
BVE coarse and are close to the DNS results, at a fraction of the computational cost. 

We start by performing a DNS computation on a fine mesh (512 x 256 spatial resolution). After a 
transient period, a statistically steady state solution is obtained at a time of around t = 10. Instantaneous 
contour plots at time t = 90 for the potential vorticity, vorticity, and streamfunction are shown in Fig. Q and 
Fig. Q for Experiment (i) and Experiment (ii), respectively. Two gyres are clearly seen in the streamfunction 
contour plots. Next, the time average of the data is taken between time t — 20 and t = 100 using 8001 
snapshots. The results are given in Fig. ([5| and Fig. We emphasize that, for both parameter sets, 
four gyres are clearly visible in the streamfunction plots. This is in stark contrast with the instantaneous 
streamfunction contour plots in Fig. Q and Fig. Q, in which only two gyres are present. The computational 
time of each DNS run was longer than 3.5 days. 
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Figure 3: Experiment (i): DNS results for Re = 200, Ro = 0.0016 and a spatial resolution of 512 x 256. Instantaneous 
field data for: (a) potential vorticity; (b) vorticity; and (c) streamfunction. Note that two gyres appear in the streamfunction 
contour plot. 
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Figure 4: Experiment (ii): DNS results for Re = 450, Ro = 0.0036 and a spatial resolution of 512 x 256. Instantaneous 
field data for: (a) potential vorticity; (b) vorticity; and (c) streamfunction. Note that two gyres appear in the streamfunction 
contour plot. 
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Figure 5: Experiment (i): DNS results for Re = 200, Ro = 0.0016 and a spatial resolution of 512 x 256. Time-averaged 
field data for: (a) potential vorticity; (b) vorticity; and (c) streamfunction. Note that four gyres appear in the streamfunction 
contour plot. 
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Figure 6: Experiment (ii): DNS results for Re = 450, Ro = 0.0036 and a spatial resolution of 512 x 256. Time-averaged 
field data for: (a) potential vorticity; (b) vorticity; and (c) streamfunction. Note that four gyres appear in the streamfunction 
contour plot. 



15 




■ ihr^i I i i i I i i i I i i i I i 
0.2 0.4 0.6 0.8 1 




0.2 0.4 0.6 0.8 1 





0.2 0.4 0.6 0.8 1 



Figure 7: Experiment (i): AD results for Re = 200, Ro = 0.0016, and a spatial resolution of 128 x 64. Time-averaged force 
function contour plots for the dissipation, subfilter-scale, Jacobian, and forcing terms. The contour intervals are between -0.05 
and 0.05 for the forcing and Jacobian terms, and -0.002 and 0.002 for the dissipation and subfilter-scale terms. 



In addition to the time series of each individual term given by Eqs. (40)-(43), following Marshall and 



Pillar 



(2011), we compute the force functions in order to analyze the contribution of individual terms, 



including the subfilter-scale term, in the new AD model (26)- (27). The defining equations for the force 



functions corresponding to the dissipation, subfilter-scale, nonlinear Jacobian, and forcing terms are: 



v 2 (/> D = A 



(45) 



V 2 (fe = S*, (46) 
V 2 0j = J, (47) 

V 2 (j) F = F. (48) 

Solving these Poisson equations for time averaged data (between t = 20 and t = 100 using 8, 001 snapshots) 
yields the corresponding force functions. Fig. ^ and Fig. ^ illustrate these force function contour plots 
for Experiment (i) and Experiment (ii) for a spatial resolution of 128 x 64. The time series of the total 
energy, dissipation, subfilter-scale, Jacobian and forcing terms are plotted in Fig. (|9|. It can be seen that 
at this spatial resolution the subfilter-scale effects are on the order of dissipation term. 
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Figure 8: Experiment (ii): AD results for Re = 450, Ro = 0.0036, and a spatial resolution of 128 x 64. Time-averaged force 
function contour plots for the dissipation, subfilter-scale, Jacobian, and forcing terms. The contour intervals are between -0.05 
and 0.05 for the forcing and Jacobian terms, and -0.002 and 0.002 for the dissipation and subfilter-scale terms. 




Figure 9: Time history of the total energy, dissipation, subfilter-scale, Jacobian, and forcing terms: (a) Experiment (i): 
S T /L = 0.04 and S M /L = 0.02; and (b) Experiment (ii): 5 T /L = 0.06 and S M /L = 0.02. 
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Next, we test the new AD model (26)-(27) using a coarser mesh with a spatial resolution of 32 x 16, 
which is 16 times coarser than the DNS mesh with a resolution of 512 x 256. On this coarse mesh, we test 
both the new AD model and the under-resolved numerical simulation BVE coarse , which does not employ 
any subfilter-scale model. To compare the models, we utilize data that is time-averaged between t = 20 and 



t = 100. The new AD model is tested with a second-order filtering operation and TV = 5 (as in Eq. (23) and 



Eq. (24)), and with the smoothing parameter a = 0.25. 

For Re = 200 and Ro = 0.0016 we plot the time-averaged streamfunction and potential vorticity contours 



in Fig. ( 10 ) and Fig. ( 11 ), respectively. We note that the new AD model yields improved results by smoothing 
out the numerical oscillations present in the under-resolved BVE coarse results. This improvement is more 



clearly displayed in the potential vorticity contour plot in Fig. (11). For Re = 450 and Ro = 0.0036 we plot 



the time-averaged streamfunction and potential vorticity contours in Fig. (12) and Fig. (13), respectively. 
The new AD model yields results that are significantly better than those corresponding to the under-resolved 
BVE coarse run. Indeed, in the streamfunction plot in Fig. ( 12 ) the new AD model clearly displays the correct 
four gyre pattern (as in the DNS plot), whereas the under-resolved BVE coarse run incorrectly yields only 



two gyres, which is nonphysical. We also plot in Fig. ( 14 ) the time history of the total energy. It is clear that 
the new AD model yields accurate results that are close to the DNS results. The under-resolved BVE coarse 
run, however, yields totally inaccurate results. The contributions of the individual terms in Experiment (i) 
and Experiment (ii), obtained using the AD model for under-resolved simulations at a spatial resolution of 



32 x 16, are also shown in Fig. (15). When we compare this to the higher resolution simulation illustrated 
in Fig. ([9]), it can be seen that the relative contribution of the subfilter-scale terms is more prominent in 
under-resolved simulations. 

To summarize, the new AD model performs significantly better than the under-resolved BVE coarse 
simulation in all tests. For the higher Reynolds number case the improvement is more dramatic. Indeed, 
the AD model correctly yields the four gyre pattern, just like the DNS computation, whereas the under- 
resolved BVE coarse computation incorrectly predicts just two gyres. We also emphasize that the new AD 
model (resolution of 32 x 16) is significantly more efficient than the DNS (resolution of 512 x 256), yielding 
a speed-up factor of more than 1, 000. 

Finally, we perform a sensitivity study with respect to the parameters used in the new AD model: M, 



the order of the high-order compact filter in Eq. (29); TV, the order of the AD procedure in Eq. (16); and a, 



the smoothing parameter of the second-order filter in Eq. (28). First, we investigate the sensitivity of the 



AD model with respect to M, the order of the high-order compact filter in Eq. (29). Since our tests show 
no significant differences among the different values of M, we are presenting below only the results obtained 



with the second-order filter in Eq. (28). 



Next, we perform a sensitivity study with respect to TV, the order of the AD procedure in Eq. (16). For 



a fixed smoothing parameter a = 0.25, the time-averaged streamfunction contour plots in Fig. (16) show a 
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Figure 10: Experiment (i): Time-averaged streamfunction data for Re = 200 and Ro = 0.0016: (a) DNS results at a resolution 
of 512 x 256; (b) under-resolved BVE coa rse results at a resolution of 32 x 16; (c) new AD model results at a resolution of 
32 x 16. Note the smoothing effect of the AD model on contour lines. Contour interval layouts are identical. 




XXX 

Figure 11: Experiment (i): Time-averaged potential vorticity data for Re = 200 and Ro = 0.0016: (a) DNS results at a 
resolution of 512 x 256; (b) under-resolved BVE coarse results at a resolution of 32 x 16; (c) new AD model results at a resolution 
of 32 x 16. Note the smoothing effect of the AD model on contour lines. Contour interval layouts are identical. 
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Figure 12: Experiment (ii): Time-averaged streamfunction data for Re = 450 and Ro = 0.0036: (a) DNS results at a resolution 
of 512 x 256; (b) under-resolved BVE coarse results at a resolution of 32 x 16; (c) new AD model results at a resolution of 
32 x 16. Note that the BVE coarse results are nonphysical, whereas the DNS and AD model results are qualitatively close. 
Contour interval layouts are identical only for (a) and (c). 




Figure 13: Experiment (ii): Time-averaged potential vorticity data for Re = 450 and Ro = 0.0036: (a) DNS results at a 
resolution of 512 x 256; (b) under-resolved BVE coarse results at a resolution of 32 x 16; (c) new AD model results at a resolution 
of 32 x 16. Note that the BVEcoarse results are nonphysical, whereas the DNS and AD model results are qualitatively similar. 
Contour interval layouts are identical only for (a) and (c). 
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Figure 14: Experiment (ii): Time history of the total energy for Re = 450 and Ro = 0.0036: DNS results at a resolution of 
512 x 256, under-resolved BVE coa rse results at a resolution of 32 x 16, new AD model results at a resolution of 32 x 16. Note 
the significant improvement of the new AD model over the results from the under-resolved BVE 

coarse run. 




Figure 15: Time history of the total energy, dissipation, subfilter-scale, Jacobian, and forcing terms: (a) Experiment (i): 
S T /L = 0.04 and 5 M /L = 0.02; and (b) Experiment (ii): S T /L = 0.06 and 5 M /L = 0.02. 
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Figure 16: Experiment (ii): Time-averaged streamfunction data for Re = 450 and Ro = 0.0036; Sensitivity with respect to 
the AD order N: (a) DNS results at 512 x 256 resolution; (b) AD model results at a resolution of 32 x 16 with TV = 5, (c) AD 
model results at a resolution of 32 x 16 with N = 3, (d) AD model results at a resolution of 32 x 16 with N = 1. Note the low 
sensitivity of the results with respect to N. Contour intervals are identical. 



low sensitivity with respect to N. The time history of the total energy plotted in Fig. (17), however, displays 



a non-negligible sensitivity with respect to N. Indeed, higher values of N yield more accurate results. Of 
course, the tradeoff is an increase in computational time due to the additional filtering operations needed for 
higher values of N: The CPU time is 170s for N = 5, 145s for N = 3, and 113s for N = 1. Thus, the two 



numerical investigations in Fig. (16) and Fig. (17) seem to suggest that, in order to balance the numerical 



accuracy and the computational efficiency, a relative low value for N might be desirable. 

Finally, we investigate the sensitivity of the AD model with respect to a, the smoothing parameter of 



the second-order filter in Eq. (28). For a fixed order of the AD procedure N = 5, the time- averaged stream- 



function contour plots in Fig. (18) show a low sensitivity with respect to a. Our numerical experiments, 
however, suggest an optimal range of 0.1 < a < 0.3, since higher values of a do not eliminate numerical 
oscillations for coarse grids. 



6. Conclusions 

This paper introduced a new approximate deconvolution (AD) model for the LES of two-dimensional 
turbulent geophysical flows. The AD model was tested in the numerical simulation of the wind-driven 
circulation in a shallow ocean basin, a standard prototype of more realistic ocean dynamics. The mathe- 
matical model employed was the barotropic vorticity equation (BVE) driven by a symmetric double-gyre 
wind forcing, which yielded a four-gyre circulation in the time mean. The AD model was tested on a mesh 
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Figure 17: Experiment (ii): Time history of the total energy for Re = 450 and Ro = 0.0036; Sensitivity with respect to 
the AD order N: DNS results at a resolution of 512 x 256 and new AD model results at a resolution of 32 x 16. Note the 
non-negligible sensitivity with respect to N. Embedded picture shows time series between t=90 and t=100. 




Figure 18: Experiment (ii): Time-averaged streamfunction data for Re = 450 and Ro = 0.0036; Sensitivity with respect to 
the smoothing filtering parameter a. (a) DNS results at a resolution of 512 x 256; (b) AD model results at a resolution of 
32 x 16 with a = 0.1, (c) AD model results at a resolution of 32 x 16 with a = 0.25, (d) AD model results at a resolution of 
32 x 16 with a = 0.45. Note the low sensitivity of the results with respect to a. Contour intervals are identical. 
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that was 16 times coarser than that used by the direct numerical simulation (DNS) run. The AD model 
yielded numerical results that were in close agreement with those of the DNS. In particular, the four gyre 
structure of the time-averaged streamfunction contour plots was recovered by the AD model. Moreover, the 
CPU time of the AD model was 1,000 times lower than that of DNS. We emphasize that this combination 
of computational efficiency and numerical accuracy achieved by the new AD model was due to the specific 
subfilter-scale model utilized. Indeed, an under-resolved numerical simulation without any subfilter-scale 
model on the same coarse mesh as that employed by the AD model produced inaccurate results. We also 
performed a numerical investigation of the sensitivity of the new AD model with respect to the parameters 
used and we found that the model is robust. This first step in the numerical assessment of the new model 
shows that AD could represent a viable tool in the LES of more realistic turbulent geophysical flows. 
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